Method of computing fir filter coefficient and program for computing same

ABSTRACT

The computer executes a first operation by a first recurrence formula, receiving a filter order (positive integer) of a universal maximally flat FIR filter, the number of zeros at z=−1 (integer equal to or more than zero), and a parameter for a group delay at z=1 (rational number). The first recurrence formula includes parameters for the filter order, the number of zeros, and the group delay, and provides coefficients in Bernstein form representation of a transfer function of a universal maximally flat FIR filter. The computer then executes a second operation composed of additions, subtractions, and division by 2 by a second recurrence formula by using a resultant of the first operation as an initial value to extract impulse response coefficients of the universal maximally flat FIR filter from a resultant of the second operation.

CROSS REFERENCE TO RELATED APPLICATION

This application is based upon and claims the benefit of priority fromthe prior U.S. provisional Patent Application No. 60/418,313, filed onOct. 15, 2002, the entire contents of which are incorporated herein byreference.

TECHNICAL FIELD

The invention relates to a computational method for impulse responsecoefficients of a universal maximally flat FIR filter, and to acomputational program for the same.

BACKGROUND ART

More than thirty years have passed since the publication of the seminalpaper of O. Herrmann on the design of maximally flat FIR digital filters[1]. Of the body of the related research that has been published sinceHerrmann's article, the most significant contribution to the theory oflowpass maximally flat FIR filters has been Baher's results on thedesign of filters with simultaneous conditions on amplitude andgroup-delay response [2]. In recent years, it has been recognized thatBaher's filters form a unifying class of maximally flat filters. Thisclass encompasses various types of seemingly distinct systems [3].Calling this class of filters the universal maximally flat digitalfilters, the authors simplified Baher's formula for the transferfunction and showed that lowpass filters of even or odd lengths, withlinear or nonlinear phase response characteristics, as well asfractional delay and half-band filters may be readily obtained bysetting the design parameters in an appropriate manner.

Not withstanding universal maximally flat digital filters have anexplicit formula for their transfer functions and, thus, may be readilydesigned using a computer algebra system, the simplified formulasprovided in [3] are still unwieldy. A direct use of the formulasinvolves computation of binomial coefficients with possibly non-integerindices that appear as terms of the summand in a three-fold sum.Binomial coefficients with non-integer indices cannot be computed by thesimple method of Pascal's triangle. Another common problem withcomputations involving the binomial coefficients is due to theirtendency to introduce very large integers to the intermediate steps ofthe computation. With the flexible computer power today, this may not bea serious impediment in many situations. However, there can be caseswhere limitations on hardware and software resources dictate the use ofefficient means in computation of filter coefficients. It is desirablethat such methods be free of binomial coefficients and have a simpleiterative structure. Moreover, it is known that if the design parameterassociated with the group delay is limited to rational numbers, as isthe case with many realistic settings, the values of the impulseresponse coefficients are rational numbers. In such cases, it isimportant to ensure that the answer to the numerical problem incomputation of an impulse response coefficient is exactly a fraction,say ⅓, not a floating number that may be shown as “0.333333335”.Therefore, it is desirable for a computation algorithm to maintain thisproperty of the impulse response coefficients and yield rational valuesfor rational design parameters.

THE FOLLOWING ARE PRIOR ART REFERENCES RELATED TO THE PRESENT INVENTION

[1] O. Herrmann, “On the approximation problem in nonrecursive digitalfilter design,” IEEE Trans. Circuit Theory, vol.CT-18, no. 3, pp.411-413, May 1971.

[2] H. Baher, “FIR digital filters with simultaneous conditions onamplitude and delay,” Electron. Lett., vol. 18, pp 296-297, 1982.

[3] S. Samadi, A. Nishihara and H. Iwakura, “Universal maximally flatlowpass FIR systems,” IEEE Trans. Signal Processing, vol. 48, pp.1956-1964, 2000.

[4] R. L. Graham, D. E. Knuth and O. Patashnik, Concrete Mathematics.Addison-Wesley, 1989.

DISCLOSURE OF INVENTION

An object of invention is to compute impulse response coefficients ofthe universal maximally flat FIR filter efficiently in a short period oftime.

According to one aspect of the present invention, the computer executesa first operation by a first recurrence formula, receiving a filterorder (positive integer) of a universal maximally flat FIR filter, thenumber of zeros at z=−1 (integer equal to or more than zero), and aparameter for a group delay at z=1 (rational number). The recurrenceformula includes parameters the filter order, the number of zeros, andthe group delay, and provides coefficients in Bernstein formrepresentation of a transfer function of a universal maximally flat FIRfilter.

The first recurrence formula, for example, is expressed asb _(j)′=(−1){(2d)b _(j−1)′+(j−1)b _(j−2)′}/(N−j+1) where 1≦j≦N withb₀′=1 and b⁻¹′=0,in which the filter order is N, the parameter for the group delay is d,the coefficients in Bernstein form representation of a transfer functionof the universal maximally flat FIR filter are b_(j)′. The resultant ofthe first operation is expressed as B′={1,b₁′, . . . ,b_(N−K)′,0, . . .,0} where the number of zeros is K.

Using a resultant of the first operation as an initial value, thecomputer executes the second operation by the second recurrence formulacomposed of additions, subtractions, and divisions by 2 to extractimpulse response coefficients of the universal maximally flat FIRfilter.

The second recurrence formula, for example, is expressed ash _(i) ^((p))=(1+E)h _(i) ^((p−1))/2+(1−E)h _(i−1) ^((p−1))/2 where1≦p≦N, 0≦i≦p with h₀ ⁽⁰⁾=B′ and h⁻¹ ^((p))={0, . . . ,0},in which a sequence for computing impulse response coefficients of theuniversal maximally flat FIR filter is expressed as h_(i)^((p))=(h_(i,j) ^((p)))=(h_(i,0) ^((p)),h_(i,1) ^((p)), . . . ), and anarbitrary sequence A_(i) is expressed as E^(j)=E (E^(j−1)A_(i)),E¹A_(i)=EA_(i)=A_(i+1), E⁰A_(i)=A_(i) in which a forward shift operatorsatisfying the expression is E. The impulse response coefficientsextracted from the resultant of the second operation are expressed ash _(i) =h _(i,0) ^((N)) where 0≦i≦N.

BRIEF DESCRIPTION OF DRAWINGS

FIG. 1 illustrates steps of introducing a computational algorithm forimpulse response coefficients according to the present invention.

FIG. 2 illustrates steps of proving the proposition 1.

FIG. 3 illustrates steps of proving the corollary 1.

FIG. 4 illustrates an overview of the computational algorithm for theimpulse response coefficients according to the present invention.

FIG. 5 is a flowchart showing the computational algorithm for theimpulse response coefficients.

FIG. 6 is a program list of the computational algorithm for the impulseresponse coefficients executed by computer.

BEST MODE FOR CARRYING OUT THE INVENTION

The embodiments of the present invention will be described withreference to the drawings in the following.

First, a description will be made on mathematical underpinning of thecomputational algorithm for the impulse response coefficients accordingto the present invention. FIG. 1 illustrates steps of introducing thecomputational algorithm for the impulse response coefficients.

It was shown in [3] that the transfer function of the universalmaximally flat digital filters, a family of filters identical to thoseproposed by Baher [2], is given by $\begin{matrix}{{{H_{N,K,d}(z)} = {\sum\limits_{0 \leq j \leq {N - K}}{{b_{j}\left( \frac{1 - z^{- 1}}{2} \right)}^{j}\left( \frac{1 - z^{- 1}}{2} \right)^{N - j}}}},} & (1)\end{matrix}$where the coefficients b_(j) are characterized by the three-termrecurrencejb _(j)+2db _(j−1)−(j−N−2)b _(j−2)=0, j≧1,   (2)with the initial valuesb ₀=1 and b ⁻¹=0.   (3)

By converting (1) to the power form representation $\begin{matrix}{{H_{N,K,d}(z)} = {\sum\limits_{0 \leq k \leq N}{h_{k}z^{- k}}}} & (4)\end{matrix}$the impulse response coefficients h_(k) become [3] $\begin{matrix}{{h_{k} = {\sum\limits_{0 \leq j \leq {N - K}}{\sum\limits_{0 \leq p \leq k}{\sum\limits_{0 \leq i \leq j}\frac{\left( {- 1} \right)^{j - i + p}\begin{pmatrix}{\frac{N}{2} - d} \\i\end{pmatrix}\begin{pmatrix}{\frac{N}{2} + d} \\{j - i}\end{pmatrix}\begin{pmatrix}j \\p\end{pmatrix}\begin{pmatrix}{N - j} \\{k - p}\end{pmatrix}}{2^{N}}}}}},{k = 0},\ldots\quad,{N.}} & (5)\end{matrix}$

Obviously, (5) is a computationally expensive formula. To amelioratethis situation, devised is an algorithm that exploits the combination of(1) and (2) to compute h_(k) for k=0, . . . ,N.

First, (1) is written in the traditional Bernstein form by introducingthe sequence b_(j)′ defined by $\begin{matrix}{b_{j}^{\prime}\overset{def}{=}\left\{ \begin{matrix}\frac{b_{j}}{\begin{pmatrix}N \\j\end{pmatrix}} & {{0 \leq j \leq N},} \\{0,} & {{otherwise}.}\end{matrix} \right.} & (6)\end{matrix}$Then obtained is $\begin{matrix}{{H_{N,K,d}(z)} = {\sum\limits_{0 \leq j \leq {N - K}}{{b_{j}^{\prime}\begin{pmatrix}N \\j\end{pmatrix}}\left( \frac{1 - z^{- 1}}{2} \right)^{j}{\left( \frac{1 + z^{- 1}}{2} \right)^{N - j}.}}}} & (7)\end{matrix}$

This is the Bernstein form of the transfer function in the traditionalsense. Although (7) has an additional binomial coefficient compared to(1), we will see later that the introduction of b_(j)′ results inconsiderable simplification.

A recurrence for b_(j)′ is obtained by substituting (6) into (2) anddividing both sides by j_(j() _(j) ^(N))≠0.

Thus, $\begin{matrix}{{{b_{j}^{\prime} + {2d\frac{\begin{pmatrix}N \\{j - 1}\end{pmatrix}}{j\begin{pmatrix}N \\j\end{pmatrix}}b_{j - 1}^{\prime}} - {\left( {j - N - 2} \right)\frac{\begin{pmatrix}N \\{j - 2}\end{pmatrix}}{j\begin{pmatrix}N \\j\end{pmatrix}}b_{j - 2}^{\prime}}} = 0},{1 \leq j \leq {N.}}} & (8)\end{matrix}$

The above recurrence can be simplified by invoking the definition of thebinomial coefficients [4] and a few simple algebraic steps. It followsthat $\begin{matrix}{{b_{j}^{\prime} = {{\frac{{- 2}d}{N - j + 1}b_{j - 1}^{\prime}} - {\frac{j - 1}{N - j + 1}b_{j - 2}^{\prime}}}},{1 \leq j \leq N},} & (9)\end{matrix}$with the initial valuesb′ ₀=1 and b′ ⁻¹=0.   (10)Thus, for given values of N and d the values of b_(j)′ can be computedrecursively by (9).

The next step is to convert (7) to the power form (4) and obtain theimpulse response coefficients. A direct expansion of the summand in (7)using the binomial theorem will introduce additional sums and newbinomial coefficients. Hence, this is not a desirable action to take.The other alternative is to use a conversion identity well known in themathematics literature and in the field of computer-aided geometricdesign. The identity is based on taking successive differences of theBernstein coefficients. We state this identity as a proposition andprovide a proof for it.

(Proposition 1)

Allow the Bernstein form of the polynomial p(x)=Σ_(0≦i≦N)p_(i)x^(i) tobe given by $\begin{matrix}{{p(x)} = {\sum\limits_{0 \leq i \leq N}{{b_{i}\begin{pmatrix}N \\i\end{pmatrix}}{{x^{i}\left( {1 - x} \right)}^{N - i}.{Then}}}}} & (11) \\{{p_{i} = {\begin{pmatrix}N \\i\end{pmatrix}\Delta^{i}b_{0}}},{i = 0},\ldots\quad,N,} & (12)\end{matrix}$where the operator Δ is the forward difference operator.

The forward difference operator is defined by the general recurrenceΔ^(j) b _(i)=Δ(Δ^(j−1) b _(i)),   (13)together withΔ¹ b _(i) =Δb _(i) =b _(i+1) −b _(i), Δ⁰ b _(i) =b _(i).   (14)

For instance, Δb₀=b₁−b₀ and Δ²b₀=b₂−2b₁+b₀. In this notation, theoperator acts on the indices, written as subscripts, of the members of asequence. To prove the proposition, we also need the forward shiftoperator E. The effect of E on the index i of a sequence like b_(i) isspecified byE ^(j) b _(i) =E(E ^(j−1) b _(i)).   (15)together withE ¹ b _(i) =Eb _(i) =b _(i+1) , E ⁰ b _(i) =b _(i).   (16)E is related to Δ byΔ=E−1.   (17)(Proof of Proposition 1)

FIG. 2 illustrates steps of proving the proposition 1. Starting with theBernstein form of p(x) and expand the term (1−x)^(N−i) using thebinomial theorem to obtain $\begin{matrix}{{p(x)} = {\sum\limits_{0 \leq i \leq N}{\sum\limits_{0 \leq j \leq {N - i}}{\left( {- 1} \right)^{N - i - j}{b_{i}\begin{pmatrix}N \\i\end{pmatrix}}\begin{pmatrix}{N - i} \\j\end{pmatrix}{x^{N - j}.}}}}} & (18)\end{matrix}$

By the symmetry of the binomial coefficients and the trinomial revisionproperty [4], obtained is $\begin{matrix}{{\begin{pmatrix}N \\i\end{pmatrix}\begin{pmatrix}{N - i} \\j\end{pmatrix}} = {\begin{pmatrix}N \\j\end{pmatrix}{\begin{pmatrix}{N - j} \\{N - i - j}\end{pmatrix}.}}} & (19)\end{matrix}$On the other handb_(i)=E^(i)b₀.   (20)

Thus, it can be written as $\begin{matrix}{{p(x)} = {\sum\limits_{0 \leq i \leq N}{\sum\limits_{0 \leq j \leq {N - i}}{\left( {- 1} \right)^{N - i - j}\begin{pmatrix}N \\j\end{pmatrix}\begin{pmatrix}{N - j} \\{N - i - j}\end{pmatrix}x^{N - j}E^{i}{b_{0}.}}}}} & (21)\end{matrix}$

The ranges for the indices of the two-fold sum above may be interchangedas $\begin{matrix}{{\sum\limits_{0 \leq i \leq N}\sum\limits_{0 \leq j \leq {N - i}}} = {\sum\limits_{0 \leq j \leq N}\sum\limits_{0 \leq i \leq j}}} & (22)\end{matrix}$with no effect on the value of the sum. Now, N−j is substituted for j inthe sum, and obtained is $\begin{matrix}\begin{matrix}{{p(x)} = {\sum\limits_{0 \leq {N - j} \leq N}{\sum\limits_{0 \leq i \leq j}{\left( {- 1} \right)^{N - i - {({N - j})}}\begin{pmatrix}N \\{N - j}\end{pmatrix}\begin{pmatrix}{N - \left( {N - j} \right)} \\{N - i - \left( {N - j} \right)}\end{pmatrix}}}}} \\{x^{N - {({N - j})}}E^{i}{b_{0}.}}\end{matrix} & (23)\end{matrix}$

The right side above simplifies to $\begin{matrix}{\sum\limits_{0 \leq j \leq N}{\sum\limits_{0 \leq i \leq j}{\left( {- 1} \right)^{j - i}\begin{pmatrix}N \\j\end{pmatrix}\begin{pmatrix}j \\{j - i}\end{pmatrix}x^{j}E^{i}{b_{0}.}}}} & (24)\end{matrix}$Note that $\begin{matrix}{{\sum\limits_{0 \leq i \leq j}{\left( {- 1} \right)^{j - i}\begin{pmatrix}j \\{j - i}\end{pmatrix}x^{j}E^{i}b_{0}}} = {\left( {E - 1} \right)^{j}{b_{0}.}}} & (25)\end{matrix}$It then follows that $\begin{matrix}{{p(x)} = {\sum\limits_{0 \leq j \leq N}{{x^{j}\begin{pmatrix}N \\j\end{pmatrix}}\left( {E - 1} \right)^{j}b_{0}}}} & (26)\end{matrix}$that is the power form representation of p(x). From the uniqueness ofthe power from representation of a polynomial and the relation betweenthe operators E and Δ it follows that p_(i)=(_(i) ^(N)) Δ^(i)b₀. (theend of the proof)

Now presented is a main and last result in this preparatory phase beforeproceeding to the derivation of the algorithm.

(Corollary 1)

The transfer function H_(N,K,d)(z) is given by $\begin{matrix}{{H_{N,K,d}(z)} = {\left( {\frac{1 + E}{2} + {\frac{1 - E}{2}z^{- 1}}} \right)^{N}{b_{0}^{\prime}.}}} & (27)\end{matrix}$(Proof of Corollary 1)

FIG. 3 illustrates steps of proving the corollary 1.

Letz ⁻¹=1−2x.   (28)Then, it can be written that $\begin{matrix}{{H_{N,K,d}(x)} = {\sum\limits_{0 \leq j \leq {N - K}}{{b_{j}^{\prime}\begin{pmatrix}N \\j\end{pmatrix}}{{x^{j}\left( {1 - x} \right)}^{N - j}.}}}} & (29)\end{matrix}$By the above proposition obtained is $\begin{matrix}{{H_{N,K,d}(x)} = {\sum\limits_{0 \leq j \leq N}{\begin{pmatrix}N \\j\end{pmatrix}\Delta^{j}b_{0}^{\prime}{x^{j}.}}}} & (30)\end{matrix}$

By substituting for x in terms of z⁻¹ in the above relation andexpanding the summand using the binomial theorem, obtained is$\begin{matrix}{{H_{N,K,d}(z)} = {\sum\limits_{0 \leq j \leq N}{\begin{pmatrix}N \\j\end{pmatrix}\Delta^{j}b_{0}^{\prime}{\sum\limits_{i \leq j}{\begin{pmatrix}j \\i\end{pmatrix}2^{- j}{\left( {- z^{- 1}} \right)^{i}.}}}}}} & (31)\end{matrix}$

By the trinomial revision identity obtained is $\begin{matrix}{{\begin{pmatrix}N \\j\end{pmatrix}\begin{pmatrix}j \\i\end{pmatrix}} = {\begin{pmatrix}N \\i\end{pmatrix}{\begin{pmatrix}{N - i} \\{j - i}\end{pmatrix}.}}} & (32)\end{matrix}$

The sum indices can also be modified to write $\begin{matrix}{{H_{N,K,d}(z)} = {\sum\limits_{0 \leq i \leq N}{\begin{pmatrix}N \\i\end{pmatrix}\left( \frac{{- \Delta}\quad z^{- 1}}{2} \right)^{i}{\sum\limits_{i \leq j}{2^{- {({j - i})}}\begin{pmatrix}{N - i} \\{j - i}\end{pmatrix}\Delta^{j - i}{b_{0}^{\prime}.}}}}}} & (33)\end{matrix}$

Using p=j−i as an index for the second sum, obtained is $\begin{matrix}{{H_{N,K,d}(z)} = {\sum\limits_{0 \leq i \leq N}{\begin{pmatrix}N \\i\end{pmatrix}\left( \frac{{- \Delta}\quad z^{- 1}}{2} \right)^{i}{\sum\limits_{0 \leq p \leq {N - i}}{\begin{pmatrix}{N - i} \\p\end{pmatrix}\left( \frac{\Delta}{2} \right)^{p}{b_{0}^{\prime}.}}}}}} & (34)\end{matrix}$

This simplifies to $\begin{matrix}{{H_{N,K,d}(z)} = {\sum\limits_{0 \leq i \leq N}{\begin{pmatrix}N \\i\end{pmatrix}\left( \frac{{- \Delta}\quad z^{- 1}}{2} \right)^{i}\left( {1 + \frac{\Delta}{2}} \right)^{N - i}{b_{0}^{\prime}.}}}} & (35)\end{matrix}$

Invoking the binomial theorem further makes it possible to write$\begin{matrix}{{H_{N,K,d}(z)} = {\left( {1 + \frac{\Delta}{2} - {\frac{\Delta}{2}z^{- 1}}} \right)^{N}{b_{0}^{\prime}.}}} & (36)\end{matrix}$By using the relation between Δ and E, (27) follows. (the end of theproof)

Next a description will be made on details of the computationalalgorithm for the impulse response coefficients according to the presentinvention.

The algorithm consists of two major stages. In the first stage, byrunning the recurrence (9) from j=1 up to j=N−K, the following sequenceis obtained.B′=(1,b′ ₁ . . . ,b′ _(N−K), 0, . . . ,0).   (37)

The entries b_(j)′ of B′ for N−K<j≦N are set to zero. The next step isto expand the transfer function using (27). This can be done by notingthat $\begin{matrix}\begin{matrix}{{\left( {\frac{1 + E}{2} + {\frac{1 - E}{2}z^{- 1}}} \right)^{p}b_{0}^{\prime}} = \left( {\frac{1 + E}{2} + {\frac{1 - E}{2}z^{- 1}}} \right)} \\{\left( {\left( {\frac{1 + E}{2} + {\frac{1 - E}{2}z^{- 1}}} \right)^{p - 1}b_{0}^{\prime}} \right)}\end{matrix} & (38)\end{matrix}$for arbitrary integer p>0. Evaluation of both sides above in the powerform yields $\begin{matrix}{{{\sum\limits_{i}{h_{i}^{(p)}z^{- i}}} = {\left( {\frac{1 + E}{2} + {\frac{1 - E}{2}z^{- 1}}} \right){\sum\limits_{i}{h_{i}^{({p - 1})}z^{- i}}}}},} & (39)\end{matrix}$where h_(i) ^((p)) is a sequence of coefficients defined byh _(i) ^((p))=(h _(i,j) ^((p)))=(h _(i,0) ^((p)) ,h _(i,1) ^((p)), . . .),   (40)and operator E effects a forward shift on the index j. Thus, it isobtained that $\begin{matrix}{{\sum\limits_{i}{h_{i}^{(p)}z^{- i}}} = {\sum\limits_{i}{\left( {{\frac{1 + E}{2}h_{i}^{({p - 1})}z^{- i}} + {\frac{1 - E}{2}h_{i}^{({p - 1})}z^{{- i} - 1}}} \right).}}} & (41)\end{matrix}$

It follows that $\begin{matrix}{h_{i}^{(p)} = {{\frac{1 + E}{2}h_{i}^{({p - 1})}} + {\frac{1 - E}{2}{h_{i - 1}^{({p - 1})}.}}}} & (42)\end{matrix}$

The ranges for i and p are specified by1≦p≦N, 0≦i≦p.   (43)

The initial values are given byh ₀ ⁽⁰⁾ =B′ ₁ , h ⁻¹ ^((p))=(0, 0, . . . ).   (44)

This completely specifies the recurrence we use in the second stage ofour algorithm. After completing the Nth step, at the point when p hasreached N, the value of the i-th impulse response coefficient is givenbyh _(i)=fist term of sequence h _(i) ^((N)) , i=0, . . . , N.   (45)

FIG. 4 illustrates an overview of the computational algorithm for theimpulse response coefficients according to the present invention. InFIG. 4, each node represents a sequence h_(i) ^((p)). The nodes markedby 0's represent zero sequences. The recurrence is illustrated in FIG. 4for N=3 and an unspecified value of K. It can be visually verified thatthe recurrence has a triangular structure.

FIG. 5 is a flowchart showing the computational algorithm for theimpulse response coefficients according to the present invention.

Receiving a filter order N (positive integer) of the universal maximallyflat FIR filter, a number of zeros K (integer equal to or more thanzero) at z=−1, and a parameter d (rational number) for a group delay atz=1, the computer computes impulse response coefficients according tothe following steps.

In step S10, all of N+1 elements B′[j] (0≦j≦N) of a single-dimensionarray B′ are initialized to zero. In step S20, all of N³ elementsr[p,i,j] (0≦p≦N, 0≦i≦N , 0≦j≦N) of a three-dimension array r areinitialized to zero. In step S30, an element B′[0] of thesingle-dimension array B′ is set to 1.

In step S40, every element of the single-dimension array B′ is decidedby changing in sequence an index j from 1 to N−K in a recurrence formulaB′[j]=(−1)×{(2d)B′[j−1]+(j−1)B′[j−2]}/(N−j+1). That is, the operation bythe recurrence formula (9) is executed.

In step S50, elements r[0,0,j] (0≦j≦N−K) of the three-dimension array rare set to elements B′[j] (0≦j≦N−K) of the single-dimension array B′.

In step S60, every element of the three-dimension array r is decided bysequentially changing, in the order of indexes j, i, p, an index j fromzero to N-p, and an index i from zero to p, an index p from 1 to N in arecurrence formular[p,i,j]=(r[p−1,i−1,j]−r[p−1,i−1,j+1])/2+(r[p−1,i,j]+r[p−1,i,j+1])/2.That is, the operation by the recurrence formula (42) is executed.

In step S70, N+1 elements h[i] (0≦i≦N) of a single-dimension array h areset to elements r[N,i,0] (0≦i≦N) of the three-dimension array r.

In step S80, the single-dimension array h is output as impulse responsecoefficients of the universal maximally flat FIR filter. This completesthe computational processing of the impulse response coefficients by thecomputer.

FIG. 6 illustrates a program list of the computational algorithm for theimpulse response coefficients by computer. A pseudo-code of a procedurethat returns the value of impulse response coefficients as an array isgiven in FIG. 6.

The following example helps elucidate the details. Let N=3, K=1 andd=(−1)/4. The actual values of the sequences h_(i) ^((p)), i=0, . . .,3, corresponding to the 3-D array r in the procedure of FIG. 6, are asfollows: h₀ ^((p)) h₁ ^((p)) h₂ ^((p)) h₃ ^((p)) p = (1, ⅙, − 11/24, 0)0: p = ( 7/12, − 7/48, − 11/48) ( 5/12, (46) 1: 5/16, − 11/48) p = (7/32, − 3/16) ( 35/48, 1/12) ( 5/96, 2: 13/48) p = ( 1/64) ( 39/64) (31/64) (− 7/64) 3:

Only non-zero values of the sequences h_(i) ^((p)), i=1, . . . ,3, arepresented. Note that all computations have been done in the field ofrational numbers in an exact manner. Hence, round-off errors areavoided.

As described above, an alternative representation for the transferfunction of universal maximally flat FIR filters enables the developmentof a very simple algorithm for evaluation of the impulse responsecoefficients of the filters. The alternative form of the transferfunction is based on the shift operator E and is derived from atechnique referred to as finite calculus in [4].

The algorithm consists of two stages. In the first stage a two-termrecurrence relation is used to obtain a 1-D sequence of numbers. If thevalue of the group delay parameter d is a rational number, then theoutputs of the recurrence are rational. The outputs are then used in asimple N-step algorithm of triangular structure. In each step of thisalgorithm, the values of the sequences obtained in the preceding stepundergo simple additions, subtractions and divisions by 2 to yieldsequences that are shorter in size. Upon the completion of the Nth step,N sequences are obtained. The first entries of these sequences are thevalues of the impulse response coefficients. As long as the computationsin the second stage are done in the field of rational numbers, the finalvalues remain rational. This provides an accurate and fast means forexact computation of the impulse response coefficients.

Here, a worst-case analysis of the computational complexity of thepresent algorithm will be described.

The computational complexity may be analyzed by evaluating thecomplexity of the simple for-loop marked by (A) and the nested for-loopdesignated by (B) in FIG. 6. For the number of additions andsubtractions C_(±), it is obtained that $\begin{matrix}{C_{\pm} = {{\sum\limits_{1 \leq j \leq {N - K}}4} + {\sum\limits_{1 \leq p \leq N}{\sum\limits_{0 \leq i \leq p}{\sum\limits_{0 \leq j \leq {N - p}}3.}}}}} & (47)\end{matrix}$For a closed-form expression, note that $\begin{matrix}{{C_{\pm} = {{4\left( {N - K} \right)} + {3{\sum\limits_{1 \leq p \leq N}{\sum\limits_{0 \leq i < {p + 1}}\left( {N - p + 1} \right)}}}}},} & (48)\end{matrix}$that can be evaluated as $\begin{matrix}{C_{\pm} = {{4\left( {N - K} \right)} + {3{\sum\limits_{0 \leq p < N}{\left( {p + 2} \right){\left( {N - p} \right).}}}}}} & (49)\end{matrix}$Thus, obtained is $\begin{matrix}{C_{\pm} = {{4\left( {N - K} \right)} + {3{\left( {{\frac{5}{6}N} + N^{2} + {\frac{1}{6}N^{3}}} \right).}}}} & (50)\end{matrix}$For the number of multiplications C* obtained is $\begin{matrix}{C_{*} = {{\sum\limits_{1 \leq j \leq {N - K}}2} = {2{\left( {N - K} \right).}}}} & (51)\end{matrix}$The number of divisions C_(/) is given by $\begin{matrix}{C_{/} = {{\sum\limits_{1 \leq j \leq {N - K}}1} + {\sum\limits_{1 \leq p \leq N}{\sum\limits_{0 \leq i \leq p}{\sum\limits_{0 \leq j \leq {N - p}}2.}}}}} & (52)\end{matrix}$that becomes $\begin{matrix}{C_{/} = {\left( {N - K} \right) + {2{\left( {{\frac{5}{6}N} + N^{2} + {\frac{1}{6}N^{3}}} \right).}}}} & (53)\end{matrix}$

In comparison, for the direct evaluation using (5), the number ofadditions and subtractions C_(±) ^(Direct) is $\begin{matrix}{C_{\pm}^{Direct} > {\sum\limits_{0 \leq j \leq {N - K}}{\sum\limits_{0 \leq p \leq k}{\sum\limits_{0 \leq i \leq j}1}}}} & (54)\end{matrix}$The right side above is of the order O(N⁴). The number ofmultiplications C*^(Direct), on the other hand, is $\begin{matrix}{C_{*}^{Direct} > {\sum\limits_{0 \leq j \leq {N - K}}{\sum\limits_{0 \leq p \leq k}{\sum\limits_{0 \leq i \leq j}{{jk}.}}}}} & (55)\end{matrix}$The right side above is of the order O(N⁵). For the number of divisionsC_(/) ^(Direct) in a direct evaluation, it can be written in$\begin{matrix}{C_{/}^{Direct} > {\sum\limits_{0 \leq j \leq {N - K}}{\sum\limits_{0 \leq p \leq k}{\sum\limits_{0 \leq i \leq j}4.}}}} & (56)\end{matrix}$The right side above, again, is of order O(N⁴). Thus, in addition to thesimplicity of most of the division operations in the present algorithm(that are divisions by 2), a significant reduction in the overallcomplexity is accommodated in terms of sheer numbers.

The invention is not limited to the above embodiments and variousmodifications may be made without departing from the spirit and scope ofthe invention. Any improvement may be made in part or all of thecomponents.

INDUSTRIAL APPLICABILITY

As described above, the impulse response coefficients of the universalmaximally flat FIR filters may be computed efficiently using a two-stagealgorithm. For a filter of order N with K zeros at z=−1, the first stageof the algorithm uses a two-term recurrence to generate a sequence oflength N−K+1. This sequence undergoes an N-step iterative process thatinvolves additions, subtractions and divisions by 2. The algorithm doesnot involve any form of binomial coefficients. Moreover, the algorithmrequires less number of arithmetical operations compared to a directevaluation using the closed-form formula. The algorithm is mostlysuitable for real-time generation of the coefficients on a DSP chipbecause it may be easily implemented in a computational environment thathas restrictions on the available hardware or software resources.

1. A method of computing FIR filter coefficients, comprising the stepsof: inputting a filter order of a universal maximally flat FIR filter, anumber of zeros at z=−1, and a parameter for a group delay at z=1, thefilter order being a positive integer, the number of zeros being aninteger equal to or more than zero, the parameter being a rationalnumber; executing a first operation by a first recurrence formula whichincludes parameters for the filter order, the number of zeros, and thegroup delay, and provides coefficients in Bernstein form representationof a transfer function of the universal maximally flat FIR filter;executing a second operation by a second recurrence formula composed ofadditions, subtractions, and divisions by 2, by using a resultant of thefirst operation as an initial value; and extracting impulse responsecoefficients of the universal maximally flat FIR filter from a resultantof the second operation.
 2. The method according to claim 1, wherein:the first recurrence formula is expressed asb _(j)′=(−1){(2d)b _(j−1)′+(j−1)b _(j−2)′}/(N−j+1) where 1≦j≦N with b₀′=1 and b ⁻¹′=0, wherein the filter order is N the parameter for thegroup delay is d, coefficients in Bernstein form representation of atransfer function of the universal maximally flat FIR filter are b_(j)′;the resultant of the first operation is expressed as B′={1,b₁′, . . .,b_(N−K)′,0, . . . ,0}, wherein the number of zeros is K; the secondrecurrence formula is expressed ash _(i) ^((p))=(1+E)h _(i) ^((p−1))/2+(1−E)h _(i−1) ^((p−1))/2 where1≦p≦N, 0≦i≦p with h ₀ ⁽⁰⁾ =B′ and h ⁻¹ ^((p))={0, . . . ,0}, wherein asequence for computing impulse response coefficients of the universalmaximally flat FIR filter is expressed as h_(i) ^((p))=(h_(i,j)^((p)))=(h_(i,0) ^((p)),h_(i,1) ^((p)), . . . ), and an arbitrarysequence A_(i) is expressed as E^(i)=E(E^(j−1)A_(i)),E¹A_(i)=EA_(i)=A_(i+1), E⁰A_(i)=A_(i) in which a forward shift operatorsatisfying the expression is E; and the impulse response coefficientsextracted from the resultant of the second operation are expressed ash_(i)=h_(i,0) ^((N)) where 0≦i≦N
 3. A program for computing FIR filtercoefficients, the program causing a computer to execute the steps of:determining every element of a single-dimension array B′ using a filterorder N being a positive integer of a universal maximally flat FIRfilter, a number of zeros K at z=−1, K being an integer equal to or morethan zero, and a parameter d for a group delay at z=1, d being arational number, all of which are provided by inputs, by changing insequence an index j from 1 to N−K in a recurrence formulaB′[j]=(−1)×{(2d)B′[j−1]+(j−1)B′[j−2]}/(N−j+1), the single-dimensionarray having N+1 elements B′[j] where 0≦j≦N, in which an element B′[0]thereof is initialized to 1 and all the elements thereof except theelement B′[0] are initialized to zero; determining every element of athree-dimension array r by sequentially changing, in the order ofindexes j, i, p, an index j from 0 to N-p, and an index i from 0 to p,an index p from 1 to N in a recurrence formular[p,i,j]=(r[p−1,i−1,j]−r[p−1,i−1,j+1])/2+(r[p−1,i,j]+r[p−1,i,j+1])/2,the three-dimension array r having N³ elements r[p,i,j] where 0≦p≦N,0≦i≦N, 0≦j≦N, in which elements r[0,0,j] thereof where 0≦j≦N−K areinitialized to elements of the single-dimension array B′[j] where0≦j≦N−K, and all the elements thereof except the elements r[0,0,j] areinitialized to zero; and extracting elements r[N,i,0] of thethree-dimension array r where 0≦i≦N as the impulse response coefficientsof the universal maximally flat FIR filter.